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Abstract 

As few real systems comprise indistinguishable units, diversity is a hallmark of nature. Di¬ 
versity among interacting units shapes properties of collective behavior such as synchronization 
and information transmission. However, the benefits of diversity on information processing at the 
edge of a phase transition, ordinarily assumed to emerge from identical elements, remain largely 
unexplored. Analyzing a general model of excitable systems with heterogeneous excitability, we 
find that diversity can greatly enhance optimal performance (by two orders of magnitude) when 
distinguishing incoming inputs. Heterogeneous systems possess a subset of specialized elements 
whose capability greatly exceeds that of the nonspecialized elements. Nonetheless, the behavior of 
the whole network can outperform all subgroups. We also find that diversity can yield multiple 
percolation, with performance optimized at tricriticality. Our results are robust in specific and 
more realistic neuronal systems comprising a combination of excitatory and inhibitory units, and 
indicate that diversity-induced amplification can be harnessed by neuronal systems for evaluating 
stimulus intensities. 
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AUTHOR SUMMARY 


Diversity is ubiquitous in natural and artificial systems, primarily arising due to spe¬ 
cialization. Specialized units have a clear role in optimally performing specific functions, 
but what is the role of the remaining non-specialized units? Surprisingly, we find that 
non-specialized units are fundamental for distinguishing the intensity of incoming inputs in 
excitable systems. Although non-specialized units themselves are inefficient at such intensity 
coding, their presence greatly enhances the performance of both the specialized units and 
the system as a whole. Our findings highlight the importance of combining both specialized 
and non-specialized units for optimal collective behavior, and indicate that diversity is more 
important than previously thought. 


I. INTRODUCTION 


[1], biological jf]] 


In numerous physical [l|, biological [2j and social [3] systems, complex phenomena (in¬ 
cluding nonlinear computations [4]) emerge from the interactions of many simple units. 
Such interactions in a network of simple (linear-saturating-response) units generate non¬ 
linear transformations that give rise to optimal intensity coding at criticality—the edge of 
a phase transition js-7]. However, optimal collective responses often require diversity js]. 
Clear examples of such optimization can be found in collective sports, business, and co¬ 
authorship in which different positions or roles require specific sets of skills contributing to 
the overall performance in their own way. 

Diversity in the nervous system, for example, appears in morphological, electrophysiologi- 
cal. and molecular properties ™ neuron types and among neurons within a single type & 
and also in the connectome [lOj, i.e., in how neurons and brain regions are connected. A large 


body of work has been devoted to show the ro_ 


topology in shaping the network dynamics 
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e of heterogeneous connectivity and network 
26]. In particular, for example, in the case 


of resonance-induced synchronization 15| . the presence or not of a single backward connec¬ 
tion may define whether synchronization or incoherent neural activity is expected in cortical 


motifs and networks, which has also been confirmed in a synfire chain configuration 27 


Crucially, diversity in the intrinsic dynamic behavior of neurons is also fundamental 


and can shape general aspects of the network dynamics 


M 


30]. The role of the inherent 
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diversity among nodes, which in many systems is at least as notable as the connectivity and 
network topology themselves, has comparatively remained largely unexplored. In particular, 
although numerous recent works have focused on optimizing features of criticality for the 
different network topologies jfi .7, 20, 21, 31-39], for convenience identical units are ordinarily 
assumed and the role of nodal intrinsic diversity on the collective behavior thus remains 
unexplored. 


Here we analyze the collective behavior at criticality in the presence of diversity in the 
excitability, which proves to be a crucial factor for the network performance. We show 
that the task of distinguishing the amount of external input, quantified by the dynamic 
range, can be substantially improved in the presence of heterogeneity. The influence of 
non-specialized units improves performance by enhancing the capabilities of both the whole 
network and of specialized subpopulations. We show the constructive effects of diversity in 
simple bimodal and uniform distributions, in more realistic gamma distributions (see Fig. [I]), 
and the robustness in networks combining excitatory and inhibitory units. 



FIG. 1. Threshold distributions in random networks. Threshold 6 indicates the minimal 
number of coincident excitatory contributions required to excite a quiescent unit. Left panel, 
bimodal distribution with 80% integrators (6 = 2). Middle panel, uniform distribution with 6 max = 
5. Right panel, gamma distribution with shape parameter a = 2, and scale parameter 6 = 1. 
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II. RESULTS 


Excitable networks with heterogeneous excitability 

Employing a general excitable model [susceptible-infected-refractory-susceptible (SIRS)], 
we characterize the dynamics and identify the constructive role of diversity in excitable net¬ 
works and neuronal systems. Node dynamics are given by cellular automata with discrete 
time and states [0 (quiescent or susceptible), 1 (active or infected), 2 (refractory)]. Syn¬ 
chronous update occurs at each time step (of 1 ms) obeying the rules: An active node j 
becomes refractory with probability 1, a refractory node becomes quiescent with probability 
7 = 0.5, and a quiescent node becomes active either by receiving external input (modeled by 
a Poisson process with rate h ), or by receiving at least contributions from active neighbors 
each transmitted with a probability A. Diversity is introduced in the threshold variable 6 J 
of each node j such that nodes with low threshold require fewer coincidental stimuli, being 
thus easily and more often excited by active neighbors than nodes with high threshold. For 
concreteness, we used Erdos-Renyi random networks with size N = 5000 and mean degree 
K = 50, but our results generalize to other sizes, connectivities, and topologies. 

Mix of specialized and nonspecialized nodes outperforms either alone 

Our analysis focuses on the input-output response function of networks subjected to 
external driving h whose intensity varies over several orders of magnitude, as is commonly 
observed in sensory systems, for example. As depicted in Fig. [2^, response functions (F) 
exhibit a sigmoidal shape with lower output rates (defined as the mean activity of the network 
or a subset thereof) for weak stimuli and high rates for strong stimuli. In this simple case 
diversity is introduced by a discrete bimodal distribution (see Fig. 1), where half the units 
are so-called integrators with 9 = 2, and the other half are nonintegrators with 9 = 1. From 
the shape of the response functions we quantify the range in which the amount of input 
can be coded by the output rate (Fig. [2^). This dynamic range A = 10 log 10 (/r 0 . 9 /ho.i) is 
a standard measure that neglects the confounding ranges of too small sensitivity [top 10% 
(F > 9 ) and bottom 10% (F < F 0 .i)], and quantifies how many decades of input h can be 

reliably coded by the output activation rate F (see caption of Fig. \2b for further details). 

Although isolated units (A = 0) code input intensity very poorly (small A), increasing 
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FIG. 2. A specialized subpopulation emerges with diversity. Bimodal distribution with 
equal numbers of integrators (9 = 2) and non-integrators (9 = 1). a, Response curves (mean 
firing rate F versus stimulus rate h) for the subpopulations of 9 = 2 (blue), 9 = 1 (red), and 
the whole network (gray). Variables Fq.i and Fq.q (red dashed lines), and ho.i and ho.g (black 
arrows) are used to calculate the dynamic range A 1 (red arrow) for the subpopulation with 9 = 1, 
where F x = F$ + xF max , h x is the corresponding input rate to the system, and Fq is the firing 
rate in the absence of input. Solid black lines correspond to the mean-field approximation (see 
Methods). Inset: Spontaneous activity Fq versus coupling strength A. b, Dynamic range A is 
optimized for different coupling strengths A for the two subpopulations. Inset: Susceptibility x 9 
for the two corresponding subpopulations; susceptibility maxima coincide with the peaks of the 
dynamic range. Susceptibility is operationally defined in the Methods. 


the contribution from neighbors (by increasing the transmission probability A) substantially 
enhances the dynamic range (Figs. Eb and [3]). However, this occurs only for coupling smaller 
than a critical value A c , at which a phase transition to self-sustained activity occurs (e.g., 
insets of Fig. Eb and Fig. |4b). As the coupling strength increases beyond the critical value, 
the dynamic range decays because the effective output range is reduced by increasing levels of 
self-sustained activity j5j. In this simple bimodal case the phase transition occurs at different 
A values for the two subpopulations, evidenced by both dynamic range A 9 and susceptibility 
X 9 (Fig. Eb and its inset). The critical value of the coupling (curve’s peak) is larger for 
integrators than for nonintegrators. Moreover, as evidenced by the difference between the 
maximum dynamic range of each subpopulation (A^^ — A aiax ~ 15 dB), nonintegrators 
greatly outperform integrators. 


5 























FIG. 3. Threshold diversity improves performance. Comparison of dynamic range for 
networks with homogeneous thresholds with 6 = 1 (green, A] 10 ™ 0 ) and 6 = 2 (purple, Al} 01110 ) with 
the 6 = 1 subpopulation of the bimodal distribution (red, A 1 ). Solid black lines correspond to the 
mean-field approximation (see Methods). 


In the presence of diversity the specialized subpopulation of nonintegrators (A 1 ) out¬ 
performs the two extreme cases with no diversity (homogeneous distribution) in which all 
units are either integrators Ai) 01110 or nonintegrators A) 10 ™ 0 (Fig. [3j). This happens because 
the response of the specialized units improves when they can also take advantage of the 
contribution of the other subpopulation of integrators, which require simultaneous neigh¬ 
boring stimulation to be effective. In the presence of integrators the network becomes more 
disconnected, requiring stronger coupling to switch to the active state. Therefore, due to a 
stronger coupling, the amplification of weak stimuli at criticality and thus the dynamic range 
are greater than in the absence of diversity. Thus, the presence of prudent units delays the 
critical transition and provides gullible units additional sensitivity to distinguish stimulus 
intensity. Remarkably, however, having all nodes behave like the specialized ones impairs 
performance. 


Tricriticality optimizes coding performance 

Henceforth, since criticality optimizes performance, we focus on characterizing the critical 
behavior for various types of diversity in the excitability. Varying the density of integrator 
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variable density of integrators 




FIG. 4. Performance is optimized at tricriticality with a critical density of integrators 
and a critical coupling strength. General bimodal distribution with varying densities of inte¬ 
grators [6 = 2) and non-integrators (9 = 1). a, Critical coupling strength (A c ) as a function of the 
density of integrators for the two subpopulations. Curves collide at a tricritical point (orange line), 
separating regimes with continuous (2nd order, green) and discontinuous (1st order, purple) phase 
transitions. Inset: Spontaneous activity Fq versus coupling strength A for the critical density of 
integrators, b, Maximum susceptibility y ma , T versus density of integrators. Inset: Susceptibility of 
subpopulation with 9 = 1 versus coupling strength for three integrator densities (0.75, 0.8, 0.85). 
c, Maximum dynamic range A max versus density of integrators. Inset: Response curves at the 
tricritical point (A = 0.1075). 


units (with 9 = 2) while the rest are nonintegrators, we find a critical point separating two 
regimes (Fig. 04): For a low density of integrators (green region) the phase transition to 
the regime of spontaneous activity is continuous (transcritical bifurcation in the mean-field 
equations for the model, see Methods); for a high density of integrators (purple region) 
the phase transition to the regime of spontaneous activity is discontinuous (saddle-node 
bifurcation in the mean-field equations) 40]. The critical coupling (A c ) grows with the den¬ 


sity of integrators for both the subpopulation of integrators (blue) and nonintegrators (red) 
and these curves collapse at the tricritical point (orange line). Apart from this collaps¬ 
ing of critical-coupling curves, the maximum susceptibility also changes qualitatively at the 
transition between the regions undergoing continuous and discontinuous phase transitions 
(Fig. 04>). Strikingly, optimal performance occurs at this transition (Fig. 04): The max- 
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imum dynamic range for generalized bimodal distributions occurs at the tricritical point 
where the sensitivity is more than two orders of magnitude larger than in the absence of 
diversity (A^ 10 ™ 0 in Fig. EJ). 


Diversity can yield multiple percolation 


Large dynamic ranges also occur at criticality in other distributions such as the uniform 
distribution. In this case, the number of units with threshold 6 is evenly distributed between 
1 and 0 mSLX , as depicted in the middle panel of Fig Q] for an exemplar case with 6 max = 
5. Notably, for the uniform distribution, A^ ax is much greater than A max of the other 
subpopulations (Fig. [5k) and of the whole network (inset). 



FIG. 5. Multiple percolation and optimal performance in uniform distributions of 
thresholds, a, Maximum dynamic range A max versus the maximum threshold of the uniform 
distribution 0 max for each subpopulation, and the whole network (inset), b, Critical coupling 
strength (A c ) as a function of 0 max for each subpopulation. The whole network (inset) exhibits 
two peaks for 0 max > 3. c, Susceptibility versus coupling strength for each subpopulation, and the 
whole network (inset). Arrows at the bottom of the panel identify the critical couplings. 


In contrast to the bimodal distribution (Fig. [Hi), the critical coupling curves of the 
subpopulations for the uniform distribution grow with 6 max without collapsing (Fig. [5h»). 
Hence, the system exhibits multiple critical couplings. However, the network taken as a 























whole exhibits only two peaks of susceptibility (insets of Figs. [5h>,c) : the lowest curve matches 
the value of the subpopulation with 6 — 1, and the other corresponds to an average of all 
subpopulations. Figure Efc displays the curves of the susceptibility for each subpopulation 
and the whole network (inset). The larger the 6 of the subpopulation, the greater the 
coupling required to optimize the susceptibility, leading to a subpopulation hierarchy. 


More realistic scenarios 

1. Gamma distribution 

The gamma distribution is more general and presumably more realistic than the bimodal 
and uniform distributions. As presented in the Methods and illustrated in Fig [U it is 
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FIG. 6. Optimal performance in gamma distributions of thresholds: the whole can 
outperform any of its parts, a, Maximum dynamic range for various subpopulations and 
the whole network (solid gray line). Inset: Gamma distribution of threshold values for shape 
parameter a = 3, and scale parameter b = 1.5. b-d, Maximum dynamic range versus scale and 
shape parameters of the gamma distribution, b, Specialized (sensitive) subnetwork; c, the whole 
network; d, difference between the whole network and the specialized subnetwork. 
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described by two independent parameters, shape a and scale factor b , and generalizes the 
exponential, chi-squared, and Erlang distributions. Exploring random networks with thresh¬ 
olds given by discrete gamma distributions (see an exemplar case in the inset of Fig. Eh), 
we find large dynamic ranges (Figs. Eh-d). The maximum dynamic range for both the sub¬ 
population with 9 = 1 and the whole network can reach ~ 40 dB. Moreover, as shown in 
Fig. Eh, the dynamic range for the whole network can outperform all subpopulations. 


2. Networks with excitatory and inhibitory nodes 


Our main result that performance can be substantially enhanced with diversity is also 
robust with respect to the presence of inhibition. Inhibition has two effects on the response 
function, influencing the dynamic range in opposite ways. On the one hand, inhibition pre¬ 
vents a rapid increase in the firing rate for small input. On the other hand, it prevents 



variable density of integrators 



density of integrators density of integrators 


FIG. 7. Robustness of optimization in a network with 20% inhibitory units. Thresh¬ 
olds are drawn from a bimodal distribution of integrators (8 = 2) and non-integrators (8 = 1). 
Maximum dynamic range versus coupling strength for the specialized subnetwork (left) and the 
whole network (right) for different types of inhibitory units: nonintegrating (triangles), integrating 
(pentagons), half integrating and half nonintegrating (squares), and the case without inhibition 
(filled circles). 
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saturation for large input. The first effect tends to reduce the dynamic range whereas the 


second effect tends to increase it. In the absence o 


a small reduction in the network dynamic range 41]. In the presence of diversity, however, 


diversity, the overall effect reported is 


we find the overall effect counterbalanced and inhibition does not alter the enhancement of 
A. Here we assume that the distribution of 6 is bimodal and 20% of the units (neurons) 
are inhibitory. After an inhibitory neuron spikes, post-synaptic quiescent neurons receive 
inhibition with probability A. Upon arrival, inhibition prevents the unit from spiking within 
a time-step period irrespective of the number of excitatory active neighbors. Figure [7] shows 
the robustness of the maximum dynamic range in the presence of inhibition. Regardless of 
whether the inhibitory units are integrators (pentagons), nonintegrators (triangles), or a mix 
of both (square) the dynamic ranges are very similar to the case without inhibition (filled cir¬ 
cles). Although inhibition has been shown to crucially shape the network dynamics 42], and 
diversity in excitatory and inhibitory populations may have different effects 43], we found 
that in the presence of diversity inhibition produces only minimal quantitative differences 
in the coding performance of networks. 


III. DISCUSSION 


Diversity has been a keystone of the recruitment theory 44] that proposed the first expla¬ 
nation for how animals can distinguish incoming input spanning many orders of magnitude, 
even when each individual sensory neuron distinguishes only a narrow dynamic range. The 
proposed mechanism there requires many neurons exhibiting responses tuned to specific 
(short) ranges of input but with the ensemble of specific ranges spanning several orders of 
magnitude. However, to satisfy this criterion sensory neurons would need to have a density of 


receptors also varying across orders of magnitude, which is not found experimentally 


44 


45]. 


A competing hypothesis claims that diversity is not required, but instead nonlinear inter¬ 


actions are sufficient 
orders of magnitude 


"or sensory systems to cope with incoming input varying over many 


46]. Remarkably, our revisited version of the recruitment theory 


reconciles the two proposals employing the key ingredient of each one: mutual (non-linear) 
interactions, which amplify the dynamic range of isolated neurons, and intrinsic diversity in 
the excitability, which requires small variability (and not variations of orders of magnitude 
as in the previous version). Therefore, showing that diversity enhances the dynamic range 


11 
















of response functions, we establish a revisited recruitment theory with solid grounds. 

Although we have focused on a specific task of distinguishing stimuli intensity, sensory 
systems also need to handle various other features. As a byproduct and another advantage of 
diversity, nonspecialized units may execute and specialize in other functions. For example, as 


recently reported in the moth olfactory system 


47], a concurrent function of the detection of 


stimulus intensity is the ability to respond promptly to external stimuli. Under evolutionary 
pressure, the ability to execute such complementary functions likely takes advantage of 
diversity to improve its own performance. 

We have demonstrated the benefits of diversity at criticality for different simple distri¬ 
butions of excitability (as requested in the recent literature 48]). Furthermore, for the first 
time we provide evidence that the well-known advantages of criticality are magnified at tri- 
criticality. The optimal performance in the simple case of two type of units is found at a 
tricritical point with a critical coupling separating the active/inactive phases and a critical 
density of integrators separating the regimes of continuous/discontinuous phase transitions. 
Even though a continuous phase transition has been proposed for the brain 7], hysteresis and 


metastability observed in models 


40 


49] and experiments 50] suggest that discontinuous 


phase transitions may also play a functional role. 

The dynamics of excitable networks exhibits two regimes: percolating (active phase) 
and non-percolating (inactive phase). As recently shown 51], percolation in core-periphery 
networks with sufficient clustering leads to double percolation, in which core nodes perco¬ 
late earlier than peripheral nodes. Analogously, for bimodal distributions we found double 
percolation with the most excitable nodes activating for weaker coupling than integrators. 
Moreover, we extended this phenomenon to arbitrarily high-order multiple percolation , with 
subpopulation thresholds following a hierarchy of excitability. 

Conclusion. Minimal models play a key role to elucidate rich emergent dynamics that 
remain elusive. Following this approach and investigating the impact of diversity in the 
intrinsic excitability, we have shown that: (i) Diversity offers clear-cut advantages in distin¬ 
guishing input with respect to homogeneous networks; (ii) At the tricritical point the system 
benefits from multiple critical instabilities, thereby optimizing performance; (iii) Subpopu¬ 
lations percolate in order of decreasing excitability; (iv) The collective response from the 
entire network can outperform all subpopulations; (v) The main results are robust to more 
realistic distributions, and can be applied to cortical systems composed of excitatory and 
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inhibitory neurons. 


IV. METHODS 

A. Network Response 

The initial condition for computing the bring rate corresponds to the active state. Nodes 
receive a strong input (h = 200 Hz) for 0.5 s, followed by a transient period of 0.5 s with the 
corresponding input level ( h ) before computing the average bring rate of each subpopulation 
over a period of 5 seconds. The reported bring rate corresponds to the average over 5 trials. 


B. Mean-Field Approximation 


In the presence of diversity the mean held map is given by a set of equations for each 


subpopulation, exhibiting a particular sensitivity to neighboring signaling 40]. For each 
subpopulation with threshold 0j, the density of refractory units R 9i at time t + 1 is given 
by R^ +1 = Ft' 1 + (1 — 7 )Rt\ where F 9, denotes the density of active units. The evolution of 
the density of active units follows F 9i = Q 9i [ 1 — (1 — h)( 1 — A/*)], where Q ( j‘ is the density 
of quiescent units, and A t 9i = Y2i= c/ (^) (AF t )*(l — A F t ) K ~ l represents the probability of not 
receiving at least 9i neighboring contributions at time f, where F t is the weighted average 
of the density d 9i of active units in each subpopulation F t = d 9i F 9i . Integrating this 
map { 52 }], we bnd the stationary distributions (F° l ) for each subpopulation. 


C. Susceptibility 

Here, susceptibility is operationally debned as x di — /(P 6i ) ~~ (P°‘)’ w h ere P° l — 

F 9i {h = 0). It was calculated over 500 trials of 100 ms after transients of 0.5 s. 


D. Gamma Distribution 

The discrete gamma distribution of thresholds is given by the smallest following integers 
drawn from the probability density function f(Q) = 6 a ~ l e~ 9 F (5“T(a)) -1 , where a and b are 
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shape and scale parameters, respectively. 
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